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Abstract 

The invaded cluster algorithm, a new method for simulating phase transitions, 
is described in detail. Theoretical, albeit nonrigorous, justification of the 
method is presented and the algorithm is applied to Potts models in two and 
three dimensions. The algorithm is shown to be useful for both first-order and 
continuous transitions and evidently provides an efficient way to distinguish 
between these possibilities. The dynamic properties of the invaded cluster 
algorithm are studied. Numerical evidence suggests that the algorithm has 
no critical slowing for Ising models. 
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I. INTRODUCTION 



This paper discusses a new cluster method for simulating both continuous and first-order 
transitions. The method, called the invaded cluster (IC) algorithm, was introduced in Ref. 
[ ij in the context of the Ising model. In this paper we give a more extended discussion of 
the IC method, present new data for Ising and Potts critical points and show how to apply 
the method to first-order transitions in Potts models. 

The Swendsen-Wang (SW) algorithms and other cluster methods |||J have led to 
vast improvements in the efficiency of simulating the critical region of a variety of spin 
models. For the Potts models, these algorithms are based on the Fortuin-Kastelyn (FK) || 
representation of the system as a correlated bond percolation problem. Each Monte Carlo 
step in a cluster algorithm consists of generating an FK bond configuration from the spin 
configuration by occupying some of the satisfied bonds (i.e. bonds across which the spins 
are in agreement) of the lattice. Clusters of connected sites are randomly and independently 
assigned a new spin value that is the same throughout the cluster. This creates the updated 
spin configuration. 

The IC algorithm shares these basic features with other cluster algorithms but differs in 
how the bond configurations are generated. For other cluster algorithms, satisfied bonds are 
independently occupied with a probability that depends on the temperature. For the invaded 
cluster algorithm, bonds are occupied in a random order until the bond configuration fulfills 
a stopping condition. For example, the stopping condition may be a requirement on the size 
of the largest cluster. For judicious choices of stopping condition the IC algorithm simulates 
the transition point of the model. 

The relation between the SW algorithm and the IC algorithm is analogous to the relation 
between ordinary percolation and invasion percolation. In ordinary (bond) percolation ||, 
bonds are independently occupied with probability p forming connected clusters of sites. At 
the percolation threshold, p c , in a large finite box, the probability that any one of these 
cluster has an extent comparable to the system size "jumps" from nearly zero to nearly one. 
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In invasion percolation 10-0], at least the version most relevant to the current work,[| bonds 
are randomly ordered and then successively occupied until a stopping condition is reached. 
If, in a system of scale L, the stopping rule is that some cluster is comparable in extent to L 
then, as is easily shown (cf. footnote (11) in [|IJ) the fraction of occupied bonds will approach 
p c as L — > oo. The IC algorithm for the g-state Potts models can be loosely described as 
the generalization of invasion percolation to the FK random cluster models. 

The IC algorithm has several very attractive features: First, the algorithm may be 
used to study a phase transition without a priori knowledge of the transition temperature. 
The IC algorithm thus enjoys the property of "uniformity" or "self-organized criticality." 
The transition temperature is an output of the IC algorithm just as p c is an output of 
invasion percolation. In cases where the critical temperature is unknown or not known with 
sufficient accuracy, this can be a significant advantage. Histogram reweighting [1^ , 14 



is 



now the method of choice for high precision measurements of the critical temperature. This 
method involves extrapolating from a guessed critical temperature and its systematic errors 
are difficult to judge. 

The IC algorithm is also an extremely fast way to simulate Ising- Potts critical points. In 
Sec. [IV D| we show that autocorrelation times for the IC algorithm are significantly smaller 
than for related cluster algorithms. Indeed, for certain quantities such as the energy and 



1 Conventionally, invasion percolation is formulated as a growth process that is initiated at a 
limited number of seed sites, e.g. a single site. 

2 In the computer science literature, an algorithm is called "uniform" if it can be applied to 
problems of arbitrary size without the need for significant precomputation. Conventional Monte 
Carlo sampling of critical points is non-uniform because precomputation is required to obtain the 
critical temperature. As the system size increases, the critical temperature must be computed to 
increasing accuracy. For the IC algorithm the critical temperature is not used as an input. In this 
setting the concept of "uniformity" is akin to the idea of "self-organized criticality." 
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the finite-size critical temperature, the integrated autocorrelation time appears to approach 
zero as the system size increases as a result of anti-correlations between successive Monte 
Carlo steps. 

Although first devised for critical points, the IC algorithm is also effective for studying 
first-order transitions. For first-order transitions, it is convenient to have the stopping rule 
control the the number of sites in the largest cluster. In this way the average magnetization is 
essentially fixed and a point in the coexistence region can be explored. By sampling a single 
point in the coexistence region the problem of exponentially long tunneling times between 



phases is avoided. The multicanonical Monte Carlo method also avoids exponential 
tunneling times and has some features in common with the IC algorithm. We also present a 
criterion for distinguishing first-order from continuous transitions. This criterion is effective 
for small system sizes and easily reveals the first-order nature of the 5-state two-dimensional 
Potts model and the 3-state three-dimensional Potts model, both of which have very weak 
first-order transitions. 

In finite volume, the IC algorithm does not sample the canonical ensemble. We call the 
invariant measure sampled by the algorithm the "IC ensemble." A fundamental supposition 
of this paper is that the IC ensemble is equivalent to the usual ensembles of statistical me- 
chanics in the thermodynamic limit. Just as the microcanonical and canonical ensembles 
agree for all local observables in the infinite volume limit, we conjecture that the IC ensemble 
agrees with these two ensembles for all local observables. On the other hand, global fluctua- 
tions differ in different ensembles. For example, energy fluctuations are proportional to the 
heat capacity in the canonical ensemble but this is not the case for the IC or microcanonical 
ensembles. Although we have no proof yet of our assertions concerning the validity of the 
algorithm, in Sec. |TTT| we carefully state our claims and supply (non-rigorous) arguments to 
back them up. 

Finally, we remark that in this study we have restricted our attention to Potts models. 
However, as shown by Wolff [Q, it is possible to generalize the cluster methods to a much 
wider class of models using an embedding procedure. Presumably, the same ideas should 



work for the IC algorithm but these matters will not be pursued here. Some additional spin 



models and graphical representations appropriate for the IC algorithm are discussed in [16 



The rest of this paper is organized as follows, in Sec. [TT] the invaded cluster algorithm 
is described in more detail and in Sec. |T| the algorithm is justified and compared to other 
simulation methods. In Sec. [TV] numerical results are presented. 



II. INVADED CLUSTER METHOD FOR POTTS MODELS 

A. Potts models 

The g-state Potts models are defined by a collection of spins, {cr^} with i belonging to 
some lattice and <7j = 1, 2 . . . q. The Hamiltonian is given by 

n = -J2(d ai>aj -i) (i) 

(id) 

where the summation is over the bonds of the lattice. If q = 2, this corresponds to an Ising 
system. Here we consider hypercubic lattices of size N = L d , usually with periodic boundary 
conditions and denotes a nearest neighbor pair. 

In the seminal work of Fortuin and Kastelyn ]5] it was shown that the above defined spin 
system gives rise to a set of percolation-type problems known as random cluster models. 
These models are defined by weights W on bond configurations u> (collections of occupied 
bonds) that are given by 

W(u) =p^(l-p)W-^q c ^ (2) 

where \u\ is the number of occupied bonds, \E\ is the number of bonds in the lattice (here 
dN), C(u) is the number of connected components of uo (counting isolated sites) and 

p = p(P) = l-e-P (3) 

is the relationship between the bond density parameter and the temperature in the spin 
system. In particular, 
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^W{u)=Z fi = Tt[e-^] (4) 

is the graphical expression for the partition function in finite volume. The weights make 
sense for any positive real q and the case q — 1 is manifestly the usual bond percolation 
problem. In this paper we only consider positive integer values of q. It has been gradually 
realized, with increasing degrees of sophistication, |p|,p!7Hl9|,H,pO[| that the graphical model 
and the spin model are "equivalent" e.g. expectations of observables in the spin model are 
easily calculable in terms of appropriate probabilities in the random cluster system. From 
the modern perspective, the two descriptions are regarded as different facets of a single 
larger problem that is none other than the annealed bond-diluted version of the original 
Potts model. 

In d > 2, for all q, there is a phase transition and for q = 1 and 2, this transition is 
continuous. (See pl |, [22] for details.) For q ^> 1, it can be proved that the transition is 
first-order [p3fl . In two dimensions, it is widely accepted that the transition is continuous 
for q = 3 and 4 and discontinuous for q > 5. In d > 3, the transition is believed to be 
discontinuous for all q > 3. 



B. The algorithm 

The invaded cluster algorithm for the g-state Potts models is defined as follows: Consider 
a finite lattice on which there is some spin configuration {cxj}. From the spin configuration, 
a bond configuration is constructed via a modified form of invasion percolation. First, the 
bonds of the lattice are assigned a random order. Bonds, are tested in this order to 

see if <7j = <jj. If the latter occurs, the bond is said to be satisfied and it is added to the 
bond configuration. These bonds are called the occupied bonds. The unsatisfied bonds are 
not considered for the remainder of the current Monte Carlo step. 

The set of occupied bonds partitions the lattice into clusters of connected sites. The 
cluster structure evolves until a stopping condition is achieved. The stopping condition is 
typically based on the size or topology of a single cluster as detailed below. When the 
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stopping condition has been satisfied, a new spin configuration is obtained by randomly 
assigning one of the q spin types to each cluster (including isolated sites) and setting each 
spin in the cluster to that value. Statistics are collected on the new spin/bond configuration 
thereby completing a single Monte Carlo step. 

Among the quantities typically measured on each Monte Carlo step is the ratio / of 
occupied to satisfied bonds. As we shall see, / serves as an estimator of the temperature of 
the system through the relation, (/) « 1 — e _/3 . 

C. Stopping rules 

The IC method does not use the temperature as an input parameter. Instead, an appro- 
priate choice of stopping rules allows us to simulate the critical point, the coexistence region 
of a discontinuous transition or a point away from the phase transition. The numerical 
results in Sec. [TV] feature the following three stopping rules: 

• Extension rule 

Some cluster has maximum extent L (the system size) in at least one of the d directions. 

• Topological rule 

Some cluster winds around the system in at least one of the d directions. (This is the 
usual rule for the identification of spanning clusters in percolation.) 

• Mass rule 

Some cluster has at least mN sites. 

It is noted that for any of the above rules, the invasion stops when exactly one cluster 
satisfies the condition. Thus, during the evolution of the bond configuration, when a given 
bond is occupied, it is only necessary to check its cluster to see if the stopping condition is 
fulfilled. The extension and topological rules involve no parameters and are used to drive the 
system to a phase transition point. Both of these rules require that some cluster is barely 
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the size of the system. We refer to these and related stopping rules as spanning rules and 
to the unique cluster which satisfies the rule as the spanning cluster. 

The mass rule involves an input parameter, m and permits us to simulate an arbitrary 
temperature and also to explore the coexistence region of discontinuous transitions. The 
mass rule is an example of a fixed parameter rule, some alternative fixed parameter rules 
that we have not yet studied numerically are: 

• Magnetization rule 

This is closely related to the mass rule and is defined in systems where the spins at the 
boundary of the sample are all set to 1. Thus bonds connecting the boundary sites to 
internal sites may only be occupied if the latter also have spin value 1. The stopping 
condition is fulfilled when the number of sites connected to the boundary first exceeds 
mN. 

• Susceptibility rule 

After the k th bond has been occupied, let us suppose that there are r = r(k) clusters, 
Ci, . . . C r containing |Ci|, . . . , \C r \ sites. The stopping condition is fulfilled when |Ci| 2 + 
. . . + \C r \ 2 first exceeds xN. 

• Energy rule 

Here one counts the number of bonds whose endpoints are in the same connected 
cluster (regardless of whether the bond itself has actually been occupied) and stops 
when the tally exceeds s\E\. 

• Density rule 

Essentially equivalent to the energy rule, we simply count the number of bonds that 
have been selected and stop when the total is p\E\. 

The magnetization, susceptibility and energy rules are derived from their thermodynamic 
namesakes via the FK representation. In a system with symmetry breaking boundary condi- 
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tions as described above, the average magnetization is precisely the number of sites connected 
to the boundary. Similarly, the susceptibility is the average size of the bond cluster contain- 
ing the origin. In periodic boundary conditions this is the same as the average of the sum 
of the squares of the cluster sizes divided by N. Finally, (cf. below) the energy per bond is 
related to the probability that both ends of a bond are in the same cluster. 

The mass rule and density rule both require some additional explanation. Let us start 
with the mass rule: in free or periodic boundary conditions, the finite-volume magnetiza- 
tion will always vanish by symmetry. Presumably, this is because all of the q distinct low 
temperature states are equally represented. Below the bulk transition temperature, this is 
manifested in the FK picture by the appearance of one large cluster and a multitude of small 
clusters. The q different values that could be assigned to the single large cluster yields the 
above mentioned convex combination. However, because of the a priori symmetry between 
these states, if we agree in advance to always assign a fixed value to the large cluster, we 
get (a finite volume approximation to) the corresponding pure state. Such a picture can be 
easily established with complete rigor if the temperature is low. For q = 2 in d = 2 this 
picture holds up to the critical point as can be derived from the main theorem in P^J23| 



and for q 3> 1 it is a consequence of the result in |2^] . Thus, we will suppose that fixing the 
fraction of sites in the largest cluster is equivalent to fixing the magnetization. 

Turning attention to the fixed bond density rule, let us assume that we are in a finite 
box A e.g. with periodic boundary conditions. The total energy divided by the number of 
bonds is given by 

£ = t4t E < 1 " <Wj >AA (5) 

1^1 <M> 

where <>/3,a denotes the thermal average at temperature in the box A. However, the 
quantity < 1 — S aha . > is easily evaluated in the FK representation as 

<l-<W rf > AA =(l-~)(l-a) (6) 

where a(q,p, A) is the probability that two neighboring sites belong to the same cluster. On 
the other hand, differentiating the expression Eq. (|2|) with respect to (3 we get 
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where b = b(p, q, A) is the probability (in a translation and reflection invariant system) that 
a given bond is present. Evidently, there is a simple relationship between the quantities a 
and 60. We thus argue that fixing the density of bonds has an equivalent effect as fixing 
the density of neighboring pairs in the same cluster. The mass and density rules have 
clear computational advantages over their counterparts, the magnetization and energy rules 
respectively. For this reason (and certain others) these rules will supersede the energy and 
magnetization rules in the remainder of this paper. 

Later we will argue that the fixed parameter rules produce, in the thermodynamic limit, 
the equilibrium system at the temperature corresponding to the chosen value of the param- 
eter. This temperature is therefore an output of the simulation. For first-order transitions, 
the goal is to find parameter values in the "forbidden region"; which is in principal easy 
for the magnetization and the susceptibility. For continuous transitions, we may probe the 
critical region by considering a sequence of magnetizations tending to zero or susceptibilities 
tending to infinity. For example, choosing \ ~ N s with s < 1 the transition temperature and 
the critical behavior of local observables will again, presumably, emerge from the output. 

The fixed parameter stopping rules are related to "constrained" random cluster ensem- 
bles. Consider, for example, the constrained mass ensemble defined by restricting attention 
exclusively to bond configurations u that satisfy the mass condition (that the largest clus- 
ter in lo has mass M = mN) and are otherwise weighted by q raised to the number of 
components of uj: 

W m (oo) = t M (uj)q C ^ (8) 

3 If 6jj is the probability that the bond < i,j > is occupied and a^j is the probability that i and 
j belong to the same cluster, the relation (q — 1)(1 — mj) = q(l — hj/p) holds even without the 
assumption of translation invariance. 
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where 1m(w) is one if u satisfies the stated mass condition and is zero otherwise. Simi- 
larly, we may construct the constrained density ensemble and the constrained susceptibility 
ensemble using weights as in Eq. [8] with 1m replaced by the appropriate restricting functions. 

Although we have not yet attempted a derivation, it would seem that standard "equiva- 
lence of ensembles" arguments could be developed to show that in the thermodynamic limit, 
these measures are identical to the usual random cluster measures at an appropriate value 
of p. In this case, according to the ideas in |L9| , the associated measure in the spin system 
converges to the corresponding Gibbs distribution. 

Somewhat to our surprise, it is readily shown that the IC density rule algorithm simu- 
lates the constrained density ensemble. This constitutes the current high water mark as far 
as rigorous results are concerned. We present the argument below: 



Theorem In finite volume, the IC algorithm with the density rule with p\E\ = P sam- 
ples the joint bond-spin distribution defined by the weights 

W(a,u) = l P (u)A((r,u). (9) 

In particular, the random cluster marginal is the constrained density ensemble defined as in 

Proof For a fixed spin configuration with at least P satisfied bonds, consider the set Qp(o~) 
of all bond configurations u consistent with the spin configuration a and with \uj\ = P: 

VL P (a) = {uj | l P (u) = 1 and A(a,u) = 1}. (10) 

We may define an "algorithm" as follows: starting from an (u),o~) with %p(u)A(a,u) = 1, 
the spin moves are identical to the usual SW or IC spin moves while the bond moves are 
defined by selecting, without bias, any u out of Qp(a). It is manifestly apparent that this 
"algorithm" samples the above joint bond-spin distribution. 

We claim that for the density rule, the bond moves are equivalent to this unbiased 
selection from Qp(a). Indeed let a denote a spin configuration. For the purpose of this 
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theorem, let us implement the IC algorithm by assigning random numbers to the bonds and 
selecting the lowest P satisfied bonds. We may write Qp(a) = . . . ,Up} with each u k 
a subset of the satisfied bonds of the a and having \io k \ = P. It is clear that under IC 
dynamics, the criterion for selecting Uk is that the highest bond in u>k is lower than the 
highest bond in any other u G Qp(a). Configurations in Qp(a) all have the same number 
of bonds, hence the probabilities are unbiased □. 



III. VALIDITY OF THE IC METHOD 
A. Comparison to the Swendsen-Wang algorithm 

The IC algorithm is in fact very similar to the SW algorithm and this, to the greatest 
extent, is the basis of our intuition concerning the former. In both cases, occupied bond 
clusters are grown on top of a spin configuration by randomly selecting some of the satisfied 
bonds. When the growth process is stopped, both algorithms generate the updated spin 
configuration from the bond configuration by the same procedure (described in the second 
paragraph of the introduction). Thus the "only" difference lies in how the bonds are selected. 

In the SW algorithm, satisfied bonds are independently accepted with probability p. If 
p is identified with a temperature as in the FK representation (Eq. |3|) one can show PJ20|| 
that detailed balance is satisfied for both the spins and the bonds: For the spins, this is with 
respect to the canonical Gibbs measure and for the bonds, with respect to the random cluster 
measure. The key observation |2(J is that the SW algorithm simulates a joint measure on 



spin and bond configurations that is defined by the weights 

W(a, oo)=p H (l- p)W-MA(a, v). (11) 

In the above, A(cr, uj) insures consistency between the bond and spin configurations: A(cr, to) 
is one if all occupied bonds in the configuration are satisfied, otherwise A(a,u) is zero. 

It is worth pointing out that the SW algorithm can in fact be described in the framework 
of an IC algorithm. Let S = S(cr) denote the number of satisfied bonds in a given spin 
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configuration. Let A p (a) denote the random variable that is chosen according to binomial 
statistics: 



The SW algorithm is defined by the stopping rule that on each round, an A p (a) is drawn, 
and then the growth stops as soon as the first A p (a) satisfied bonds are accepted. Viewed 
in this light, the SW algorithm is seen to be quite similar to some of the IC algorithms 
under consideration. However, unlike the SW algorithm, we are unable to write down a 
joint distribution similiar to Eq. |TT] for any of the IC ensembles except for the case of the 
density rule. 



As the title of this subsection indicates, we will provide separate discussions of the critical 
and non-critical algorithms. In part, this is because the latter require fewer assumptions 
about the equilibrium state of the Potts/FK system. For the fixed parameter stopping 
rules, the fundamental assumption is that the ensemble sampled by the IC algorithm has 
the same local properties as the canonical ensemble at a temperature that yields the chosen 
value for the thermodynamic function (e.g. magnetization or susceptibility) that defines the 
stopping rule. We present the argument for the susceptibility rule though similar arguments 
are possible for other rules. 

Consider a large Potts/FK system in equilibrium at p < p{(3 c ) (i.e. T > T c ). The "system 
wide" average cluster sizes will be very close to the mean value which, in turn, is close to 
the limiting infinite volume susceptibility, xip)- Explicitly, for bond configuration ojl on a 
lattice of scale L, we may compute 



P(A p (a)=A) = ( s A )p A (l-p) 



S-A 



(12) 



B. Away from criticality: fixed parameter algorithms 



CVl) 



CI + ... + c f 2 



(13) 



N 



and, as L — > oo, the distribution of C(u>l) will be sharply peaked about xip)- 
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Now let us contrast the behavior of the SW and IC algorithms on a given spin configu- 
ration. Suppose that both the SW and IC algorithms run by assigning a random number 
uniformly in [0, 1] to each bond of the lattice. For the IC algorithm, this provides us with 
the ordering while for the SW algorithm at parameter p, the instructions are to occupy all 
satisfied bonds with random number less than p. For the SW algorithm we may occupy 
bonds as a function of continuous time t, < t < 1 . At time t, all satisfied bonds with 
value less than t are occupied. Clearly, if we were to stop too soon, e.g. at t = p — e, the value 
of C will be, with high probability, strictly smaller than the equilibrium value. Similarly if 
we were to go beyond p to t = p + e, we will get a value of C that is too large. On the other 
hand, stopping, at t = p as we are supposed to, yields a value of C that, by definition, is 
typical of the equilibrium distribution. 

Now, starting from the same spin configuration, let us do a step of the IC algorithm 
(with the susceptibility rule) stopping when C = x — x(p)- We may also envision this 
operation taking place as a function of continuous time. Now at time t, the same bonds 
have been collected in both Monte Carlo schemes and we can reiterate the previous discussion 
to conclude that the algorithm will not stop significantly before or after t = p. If the system 
is in equilibrium initially then the IC algorithm chooses nearly the same stopping time as 
the SW algorithm. Since the latter was in equilibrium, evidently under IC, we stay in 
equilibrium. 

The opposite perspective provides us with an equally valid argument: Suppose it is the 
case (as is observed) that the fraction / of satisfied bonds that are occupied has a distribution 
that is sharply peaked. Then, each step of the IC algorithm amounts to an iteration of the 
SW algorithm with parameter value equal to /. If the stopping rule demands that C = xip)i 
it follows that / = p. 

Of course for an actual simulation in finite volume, the above arguments are by no means 
a rigorous proof: To achieve C = x{p)-> the algorithm stops at t — p + tjl where i]l is a 
random variable depending on the spin configuration. Even if we know that 1 77^, j — > as 
L — > 00, we cannot, as of yet, control the effects that these fluctuations have on the limiting 
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distribution for the IC algorithm. Indeed, the IC distribution will differ from the canonical 
distribution in finite volume. However, it is hard to believe that these objects do not tend 
to the same distribution in the infinite volume limit, all we lack is a proof. Nevertheless, to 
summarize our argument, the assumption that the bond/spin configurations typical of the 
Gibbs distribution are close to the ones of the IC ensemble is self-consistent in and of the 
fact that iterations of the IC algorithm keep us in the vicinity of the Gibbs distribution. 

Let us now turn our attention to a discussion of the situation when the spin configuration 
is not typical of the desired equilibrium distribution; here the arguments will be somewhat 
less complete. We will again consider a given spin configuration and compare what happens 
under SW vs. IC dynamics. Suppose, for example, that the spin configuration is at too 
high a temperature. We may imagine that we have a configuration that is typical of the 
Gibbs distribution at an inverse temperature (3 such that 1 — e _/3 = p < p = 1 — e _/3 . If 
we do a single iteration of the SW algorithm, again collecting our bonds according to the 
time parameter t, as always, we stop when t = p. Because the temperature of the spin 
configuration was higher than 1//?, the percentage of satisfied bonds will be relatively lower 
then typical for configurations that are at the right temperature. Evidently, when we stop, 
with large probability the average bond cluster size will be smaller than xip) but larger than 
x(P)- The resulting spin configuration will therefore be of intermediate character between 
those that are typical of temperatures 1/(3 and temperatures 1/(3. By contrast, when t — p, 
the invaded cluster algorithm does not stop. Hence, the fraction / of satisfied bonds that get 
occupied under IC dynamics is in excess of p and the new configuration (also of intermediate 
character) is further towards the low temperature side. It is thus apparent that away from 
equilibrium, we are pushed in the right direction and, as this example illustrates, we are 
pushed harder in the direction of equilibrium under IC dynamics than under SW dynamics. 
Similar arguments apply to the case p < p and show that now / < p. We believe, this 
"negative feedback mechanism" is ultimately responsible for the immense reduction (or 
complete absence) of critical slowing down that results from using IC dynamics. 

Unfortunately, we now leave terra firma in order to discuss the case where the configu- 
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rations cannot be characterized in terms of a single temperature-like parameter. In general, 
dynamically generated configurations are out of equilibrium, however, we can consider the 
following proposal for the definition of an effective temperature: For a given configuration 
<jl let C(aL,t) be the mean square cluster size (as defined in Eq. |13|) observed when we 
have occupied all satisfied bonds with value less than t. Let C(aL,t) be the average of 
this quantity over all realizations of random numbers on the bonds. A measure of the ef- 
fective temperature of the configuration is to compare C(ct.l,£) with xip)i the actual 
(equilibrium) susceptibility as a function of the temperature parameter p. If there is a single 
(non-zero) point where these curves cross, this may be identified as an effective temperature. 
In any case, if C(<jl, t = p) is smaller than x(p) we can assert that the temperature is higher 
than that corresponding to the parameter p while if C(ai,t = p) > x(p) ^ is lower. Ob- 
viously the colder configurations should be heated up and the hotter configurations should 
be cooled down. For spin configuration a and target susceptibility Xi the typical simulation 
temperature of a move will be given by t(x) such that C(aL,t(x)) — xip)- Since the curve 
C(o"l, t) is monotone we again see the negative feedback mechanism-if the spin configuration 
is colder than that corresponding to p the simulation temperature will be higher than that 
corresponding to p and vice versa. 

Of course these considerations apply as well to the other fixed parameter IC algorithms 
where we would argue-just as persuasively-that the magnetization vs. t or energy vs. t 
profiles of a spin configuration can be used to measure an effective temperature. It is worth 
remarking that there is one system where these suppositions are exact, namely the long- 
ranged mean-field Ising model. Here, for an N site system, each site interacts with every 
other site via a coupling that scales inversely with N. At non-critical temperatures, the SW 
algorithm drives this system to equilibrium exponentially quickly (as expected) while, by 
contrast, the IC algorithm achieves equilibrium in at most two Monte Carlo steps [|1(| . 
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C. At criticality: spanning algorithms 



The general philosophy that underpins our belief in the validity of the critical algorithms 
is quite similar to the non-critical cases. The important differences lie in the implicit assump- 
tions we have made concerning the behavior of the graphical representation at criticality -in 
particular in finite volume. Indeed, for a Potts system with a continuous transition, if we 
ask for the value of the parameter where the susceptibility is equal to a certain value, the 
answer is unambiguous. Furthermore, provided that L is large compared to the typical 
size of bond clusters, the statistics in a finite lattice of side L should represent an excellent 
approximation to the infinite volume behavior. 

For the critical algorithms, the entire premise begins with non-trivial questions about 
the equilibrium critical behavior of the random cluster model in finite volume. For example, 
in a large system is there a single cluster with the following two properties? 

1. The extent of the cluster is the scale of the system. 

2. The cluster does not contain two (or more) disjoint subsets each of which satisfy 
condition (1). 

Under the standard assumptions concerning the nature of the critical point [] the following 
picture, in the graphical representation, emerges: Above T c , the probability that there is 
any such cluster goes to zero exponentially in the scale of the system. Below T c , there is a 
single large cluster that exhausts a fraction-equal to the spontaneous magnetization-of the 
system. Within this cluster, there are many separate paths that are the scale of the system. 

4 These include absence of an intermediate phase, uniqueness of the phase above and at T c , unique- 
ness of the q magnetized phases among the translationahy invariant states below T c , continuously 
diverging correlation lengths using several definitions of the correlation length and positivity of 
the surface tension below T c . Nearly all of these properties have been established for independent 
percolation and a significant (but not sufficient) fraction in the Ising case. 
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To prevent this large cluster from happening, requires a fluctuation presumably as rare as 
exp[— const. L d_1 ]. By process of elimination, the only place that a cluster with the above 
properties could exist on all scales is the critical point. 

This is not to say that the above line of reasoning proves that these or other kinds of 
spanning clusters are indeed typically observed at the critical point (e.g. they could be power 
law rare) however, as we will see, the full validity of such an assertion is not essential for 
the broader features of our argument. 

Indeed, the basic reasoning is now pretty much the same as in the non-critical case: If 
the distribution of / (the fraction of satisfied bonds that are occupied) is sharply peaked, the 
central value must correspond to that of the critical parameter, p(/3 c ). If this central value 
were too small then the correlation length would be a small fraction of the system and the 
clusters would not get big enough. If the peak value is so large that the (low temperature) 
correlation length is small compared to the system size, the biggest cluster would be too 
big. Thus, the value of / has to be close enough to the critical point to ensure that the 
correlation length is at least a scaling fraction of the system size. 

If the required spanning cluster is itself, somehow, atypical of criticality, the above ar- 
gument is still valid. As a concrete example, consider a stopping rule that terminates the 
cluster growth when there is a cluster of size VT. Such a bad choice of a spanning cluster 
will nonetheless heavily favor the critical value of / over any non-critical value. The worst 
that could happen is that the scaling of the spanning cluster itself might be of the wrong 
type but in any case all local observables will still take on their critical values. Finally, 
starting at non-critical spin configurations, the negative feedback mechanism discussed in 
the previous subsection applies to these algorithms as well. 

The weaker point in our argument concerns our reasoning as to why the distribution 
of the / values should be sharply peaked. (First and foremost, this has been observed 
in every critical system.) Let us imagine the problem in an infinite volume setting and 
consider a critical spin configuration. We again regard the process of growing the clusters 
as an independent percolation problem defined on the random graph that is provided by 
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the satisfied bonds of the spin configuration. Let us first assume that, in the usual sense, 
this problem has a sharp percolation threshold, t c . It then follows easily that t c corresponds 
to the critical value of the FK parameter; t c = 1 — e^ c = p{(3 c ). Indeed, at t = p(/3 c ), we 
have achieved the critical FK bond configuration and our assumption of a sharp threshold 
rules out the possibility of any other value. Going back to finite volume and starting from a 
critical configuration, the argument in this case is finished: If we stop at an / significantly 
different from p(/3 c ), we will get the wrong sort of clusters and stopping at / ~ p(/3 c ) we 
keep the spin configuration critical. 

Unfortunately, in d = 2, it is not the case that the underlying percolation problem has 
a sharp transition. Specifically, in the Ising model on the two-dimensional square lattice 
it was shown in [^7],^!| that percolation of one spin type is necessary and sufficient for the 
existence of a positively magnetized phase. (The analog of this result for the general g-state 



Potts models was proved in p9| ) 

In particular, this means that in a critical configuration, there is no infinite cluster of 
satisfied bonds and thus, even if t = 1, there is no percolation in our secondary process. 
Evidently the percolation clusters on the critical spin configuration will themselves look 
critical for all t between p(T c ) and 1. It is easy to believe (but hard to prove so we will 
spare the reader the details of the heuristics) that the clusters will not go critical until 
t = p(T c ). Thus, the algorithm will not stop collecting bonds until at least this point. 
However, in finite volume, there are, undoubtedly, typical critical spin configurations that 
forbid the existence of a spanning cluster. Q Of course this kind of disaster is ruled out by 
the mechanics of the algorithm: whatever the stopping condition, if it was satisfied on the 



5 For example, if there are star-connected chains (meaning that neighbors and next nearest neigh- 
bors both count as connected) of plus spins and of minus spins winding both ways around the torus, 
the topological condition cannot be satisfied. At the critical point, such configurations presumably 
have uniformly positive probability on all scales. 
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last iteration, it is satisfiable on the next one. However, near disasters can occur causing 
"bottlenecks" -situations where one of a relatively few bonds must be occupied in order to 
achieve a spanning cluster. This would have a tendency to drive us to higher values of /. 

We believe that these bottlenecks do occur and, in fact, are responsible for the relatively 
broad tails in the distribution of / in the region p(f3 c ) < f < 1 that have been observed 
in our two-dimensional simulations. However, we also believe that these events affect only 
the details of how the L — > oo limit is achieved, not the limit itself since there are alternate 
routes circumventing bottlenecks occurring on all scales. Nevertheless, the finite-size scaling 
is sometimes quite complicated and, in certain instances, we must resort to semi-empirical 
fitting of the data. 

An interesting feature of the IC algorithm (using spanning rules) is that the approach 
to equilibrium is along a critical trajectory. For example, if the starting configuration is 
characteristic of zero temperature, the initial bond configuration is typical of ordinary bond 
percolation at threshold. Thus some sort of power law correlations are actually established 
on the first step. 

IV. NUMERICAL METHODS AND RESULTS 

A. Implementation of the algorithm 

The most difficult part of the IC algorithm is the construction of a cluster configuration 
from a spin configuration. The first step is to produce a random permutation of the set, E 
of bonds of the lattice (\E\ = dL d here). This is accomplished through \E\ random pairwise 
permutations. Initially, let tx : {1, . . . , \E\} — > E be some conventional initial order on E. 
For j — 1 to l-El, 7T is updated by choosing a random number, r in the range j to \E\; 
then the j th and r th elements of the permutation are interchanged, ir(j) # ir(r). It is well 
known that after \E\ steps, 7r is a random permutation. The computational work involved 
in making the random permutation is nearly linear in the number of bonds. 
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Bonds are explored in the order given by the random permutation. If a bond is satisfied 
it is added to the cluster configuration. The data structure for the cluster configuration and 
its updating is done in the same general way as for other cluster and percolation algorithms 
using the Hoshen-Kopelman (or "disjoint set forests") method. Each cluster of sites is 
described as a rooted tree and when two clusters are combined, the root of the smaller 
cluster becomes a son of the root of the larger cluster. When two clusters of the same size 
are combined the conventional direction associated with the bond that joined the clusters 
determines which site is to be the root. Information concerning the current state of the 
cluster as a whole, such as its mass, is stored with the root. 

After the cluster configuration is updated by the addition of a bond it is necessary to 
check whether the current cluster satisfies the stopping rule. For the fixed parameter rules 
this is straightforward. For the other spanning rules, it is necessary to associate with each 
site a vector from the site to the root of its cluster. The set of distance vectors, {vi} is 
updated in the natural way when two cluster are combined. The sites in the larger cluster 
retain their previous coordinates relative to the root. The sites in the smaller cluster take 
new coordinates, v', 

v'i^-vi- v k - e jk + vj (14) 

where ejk is the unit vector of the new bond added to the lattice, j is the site in the larger 
cluster which connects to k in the smaller cluster. For the topological rule, stopping can 
only occur if the new bond is added as an internal bond. If the new bond (j, k) is an internal 
bond we evaluate, 

v l «- Vj - e jk . (15) 

If vl 7^ v k the current cluster is multiply connected and the topological rule is satisfied. 

For the extension rule, each cluster must have associated with it the coordinates of the 
2d sites which are the most distant from root along the d axes in the positive and negative 
directions. Updating these coordinates after two clusters are combined is somewhat tedious 



due to periodic boundary conditions and is described in |30 . 
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The invaded cluster algorithm requires 1.8xl0~ 5 seconds per update per spin on a DEC 
Alpha 2100 workstation. The running speed is about a factor of two slower than for the SW 
algorithm. 

For the results reported here we start with an initial ordered (T = 0) configuration. 
Unless otherwise stated, the system is allowed to equilibrate for two hundred Monte Carlo 
steps (MCs) before data collection. If no error bars are shown in a figure, the error is smaller 
than the symbol size. 



B. Continuous transitions 



1. Three-dimensional Ising model 



Figure [I] shows data for the mean value of the ratio of occupied to satisfied bonds, (/) 
for the three-dimensional Ising model as a function of L -1 ' 59 . The power of L is chosen to 
approximate the inverse of the three-dimensional Ising correlation length exponent \jv « 
1.59. Results for both the topological and extension spanning rules are shown. Finite-size 
corrections are smaller for the topological stopping rule. The best linear fit to the data for 
the topological rule yields 0.35803 in comparison with a recent value p(T c ) =0.358098(7) 



from Ref. 0. 



Using ideas made plausible by finite-size scaling theory we can obtain two independent 
critical exponents. Figure || shows log((/) 4 — (f) e ) plotted against log(L) for the three- 
dimensional Ising model where (f)t is measured using the topological rule and {f) e is mea- 
sured using the extension rule. A fit to the data yields v = .63, which is in agreement with 



the value 0.6289(8) from A second independent exponent can be obtained either from 
the cluster size distribution of the scaling of the largest cluster. Figure |] shows n(s), the 
number of clusters per site of size s with the data binned in octaves. This figure shows that 
the cluster size distribution is indeed self-similar and allows us to estimate the exponent r, 
defined by n(s) ~ s~ T . A straight line fit yields r = 2.19 compared to the accepted value 
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of 2.21. A more efficient way to obtain a second independent exponent is via the fractal 
dimension of the spanning cluster. The average size of the spanning cluster is plotted against 
the system size in the inset of Fig. from which we obtain (3/v = 2.45 compared to the 
accepted value 2.47. 

2. Critical two-dimensional Potts models 

Figure ^ shows results for the extension rule applied to the two-dimensional Ising model. 
Both the mean and median value of / are plotted against L~ l (in accord with finite-size 
scaling since v — 1 for the two-dimensional Ising model) and are seen to converge to the 
exactly known value of p(T c ) . The fact that median lies below the mean shows that the dis- 
tribution is skewed toward larger values of /. This is presumably caused by the simultaneous 
percolation of the spins discussed in Sec. [Ill (J| . This is in contrast to the three-dimensional 



Ising model for which the / distribution is very symmetrical. The inset to Figure |] shows 
var(Z) 1 / 2 plotted against 1/L. The solid line is a fit to the data whose leading behavior 
is L -1 / 2 . This curve supports the hypothesis that the distribution of / becomes sharp as 
L — > oo. 

Figure ||] is a plot of the average energy per spinF] vs. L" 1 with the exact value plotted on 
the vertical axis. A fit to the data of the form, Eq + EiL" 1 + E2L~ 2 yields, Eq = —1.706 which 
is reasonably close to the exact value —1.7071 .... Energy fluctuations are shown in the inset 
of Fig. The quantity var(e)iV is seen to increase roughly linearly in L. This is in contrast 
to the canonical ensemble where var(e)iV is the specific heat and diverges logarithmically in 
L for the two-dimensional Ising model. The behavior of energy fluctuations underscores the 
difference between the IC ensemble and the canonical ensemble. 

The top panel of Fig. |6| shows the fraction of occupied bonds vs. the mass of the largest 
cluster for the mass rule. Data collapse for a range of system sizes predicted by the finite-size 



3 In this section, e refers to the energy per spin rather than the energy per bond. 
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scaling ansatz, 

[{f)-p(T c )]L^~G(mL^). (16) 

is confirmed in the lower panel. These results demonstrate that the IC algorithm can be 
used to extract quantitative results for the critical temperatures and critical exponents using 
systems of modest size. 

Figure [7] is a log-log plot of the deviation of / from its exact value versus the system 
size for two-dimensional Potts models with continuous transitions, q = 1,2,3,4. Except for 
the Ising case we have used the topological rule. The extension rule is used for the Ising 
case. Figure H is a log-log plot of var(Z) 1 / 2 vs. the system size for the two-dimensional 
Potts models with continuous transitions. The figure shows that the / distribution becomes 
narrow as a power of the system size L. Fitting the last five data points for each q to the 
form, var(Z) 1 / 2 ~ L~ b ^ yields 6(1) = .71(1), 6(2) = .46(2), 6(3) = .30(2) and 6(4) = .23(1). 
For percolation this result is close to the expected scaling, 6(1) = l/u(l). For invasion 
percolation it is believed that the full / distribution scales with L~ l / U ^\ For the other values 
of q, 6(g) is much smaller than l/v(q) and decreasing with q. We do not yet understand the 
finite-size scaling of the / distribution. 

C. First order transitions in Potts models 

1. Three-dimensional 3-state Potts model 

We first discuss results for the topological and extension stopping rules for the three- 
dimensional 3-state Potts model. This model has a weak first order transition. Fig. |9| 
shows the mean and median values of / vs. L^ 1 for system sizes between 10 and 70 for 
the topological stopping rule and between 10 and 110 for the extension stopping rule. Data 
was taken from samples of 10 MCs for each lattice size up to L = 70 and between 3000 
and 6000 MCs for the larger sizes. In the case of the topological stopping rule, a fit of 
the mean to a function of the form cq + C\L~ X I 2 + C2L' 1 + c^L' 2 gives Cq = 0.4232. This 
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is in good agreement with results obtained with other methods [j32 |. With the topological 
stopping rule, the median shows no finite-size effects within the error, as is the case for 
the three-dimensional Ising model. It can be fitted to a practically horizontal line, which 
extrapolates to 0.4228. Again the median seems to be the better choice if predictions about 
the infinite system are to be made, although it tends to be more noisy. As in the case of the 
two-dimensional Ising model, the difference between the mean and the median results from 
a tail in the distribution of / towards / = 1. 

For the extension stopping rule, the median enters a flat region starting at about L = 40. 



The arithmetic mean of the last seven data points is .42336. Current values of T c |32| agree 
in the first four digits, and so does the value that we obtained from the median in this way. 
Again, the mean is above the median and starts to approach it only for very large system 
sizes (L > 80). 



2. Mass rule for first-order transitions 

The above results show that spanning rules may be used to accurately locate a weak first- 
order transition however, they perform poorly for strong first-order transitions (large values 
of q). The difficulty is that the / distribution becomes increasingly broad and asymmetric 
with a tail extending toward / = 1. We believe that the tail in the / distribution is related 
to the way the spanning condition is met for strong first-order transitions. The spanning 
cluster is a nearly linear object which extends across the system in a background of small 
clusters whose size is presumably the correlation length. For large values of q the spanning 
cluster is very narrow (somewhat like a river running through a terrain of small clusters). 
This observation is consistent with the increase in p(T c ) with increasing q. In addition, for 
large q, there are severe bottlenecks; one of only a few bonds must be occupied to meet the 
spanning condition. This leads to a broad / distribution with a tail toward / = 1. Although 
we believe the / distribution becomes sharp for large L the convergence is very slow. An 
additional difficulty in using either of the spanning rules for strong first-order transitions 
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arises from the fact that the spanning cluster is nearly reproduced in successive Monte Carlo 
steps so that the autocorrelation time is large. 

These problems can be avoided with the mass stopping rule. At the phase coexistence 
temperature T c , the magnetization may take any value from to mi where mi is the magne- 
tization of the pure low temperature phase. Thus we expect / to approach p(T c ) for every 
m between and mi. At mi we expect the derivative of f(m) to jump to a finite value as 
the systems enters the ordered phase with T < T c . 



Figure [10] shows a plot (/) vs. m for the two-dimensional 10-state Potts model, obtained 
from the mass rule for system sizes 50, 100 and 200. Data was obtained after an equilibration 
of 1000 MCs from a sample of 2 x 10 4 to 5 x 10 4 MCs. The dashed line denotes the exact value 
of p(T c ) = y/q/(l + y/q). It is clear that the crossing point of these curves (near m = 0.6) 
for different system sizes provides an accurate estimate of p{T c ). Furthermore, the curves 
become increasing flat for large L and presumably converge (non- uniformly) to p(T c ) . The 
value of mi can be estimated from the largest value of m for which (/) = p(T c ). From data 
for L = 50, 100,150, 200 and 500 (the data for 150 and 500 is omitted from the plot for 
clarity) we find convergence to the value, mi = 0.8544. 

For small values of m we find a region which moves increasingly toward zero where / is 
significantly greater than p(T c ). The non-monotonicity of / as a function of m occurs only 
for those Potts models with first-order transition (compare Fig. |i~2"l) . 

In the inset of Fig. [H^ the standard deviation of / is plotted vs. 1/L. This quantity 
appears to vanish as L — > oo though at a rate that depends on m. The data is consistent 
with our belief that the / distribution approaches a delta function at p(T c ) for any m < m\. 
Any mixture of low and high temperature phases in the coexistence region can be sampled by 
fixing the ratio, m/mj. This argument can be confirmed by looking at the energy per spin. 
Let ni be the fraction of the system that is in the low temperature phase and assume that 
it is proportional to m. If ei (e^) is the energy per spin of the pure low (high) temperature 
phase, we expect the energy per spin e to behave like 
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(s) = {m/mi)ei + (1 - m/ mi )e h + e sf A/N (17) 

at the transition point. The third term includes the interfacial energy per unit area e s f and 
the area A and should vanish with system size like 1/L. 

Fig. |TT| shows (s) vs. m for L = 50, 100 and 200 from the same runs as the data of Fig. 



TU| . Error bars obtained with the jack-knife method are smaller than the symbols. There is a 
large region from small m to about m\ = 0.855 where the energy is described by a line with 
negative slope plus a correction that vanishes as L becomes large. Again, this statement 
is also based on additional data for L = 150 and 500. This behavior is in good agreement 
with Eq. [17|. At mi the energies for different L collapse, as there is no interfacial energy left, 
and the systems enter the ordered phase below T c , causing the energy to drop. Only for 
small m the energies leave the presumed curve due to finite size effects. In order to estimate 
the energy of the high temperature phase, we fitted the L = 200 data points in the range 
0.05 < m < 0.5 to a function of the form cq + C\m + c^vn 2 , that gives Eh = Cq = —0.969, 
which agrees well with the exact result —0.9682 . . . P3"| . Evaluating the data for the biggest 
system used, L = 500 at mi gives e\ = -1.661, which is close to the exact value —1.6642 . . . 
as well. 

There are few computational methods for reliably distinguishing the order of a phase 



transition f34jl . The 3-state Potts model in three dimensions for example at the phase 
transition point in many respects behaves just like a second-order transition. If the magne- 
tization rule is used, however, even weak first-order transitions seem to behave differently 
than continuous ones. We observed that / is a non-monotonic function of m for first-order 



transitions. Figure [12] shows the median value of / vs. m for Potts models for several values 
of q. The known values of p(T c ) are marked by dashed lines. Note that the curves are mono- 
tone increasing for models with second-order transitions, (q — 2, 4) and non-monotonic for 
models with first-order transitions (q — 5, 6 and 10). Non-monotonicity is also found for 
the three-dimensional 3-state Potts model, see Fig. |TB|. It should be noted that both the 
three-dimensional 3-state and two-dimensional 5-state Potts models have extremely weak 
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first-order transitions so that this criterion is quite sensitive. It is useful even if the correla- 
tion length is larger than the system size. 



D. Dynamics of the IC algorithm 

In this section we study the dynamic properties of the IC algorithm. The normalized 
autocorrelation function of an observable A is defined by, 

„ . , < A A t > - < A > 2 . . 

< A 2 > - < A > 2 



where t is time in Monte Carlo steps. The three sets of points in Fig. 13 are the normalized 
autocorrelation functions of the absolute value of the magnetization m, energy e, and fraction 
of occupied bonds / for the two-dimensional Ising model. Numerical data were collected for 
the topological stopping rule from a run of 10 4 MCs which was divided into 10 groups with 
errors estimated by the jack-knife method. In a few steps all three autocorrelation functions 
have nearly vanished. The autocorrelation functions of / and e display a negative overshoot 
on the first step which becomes larger for larger system sizes. 

The integrated autocorrelation time ta, which is required for estimating the errors in 
measuring the observable A, is defined by, 

TA(w) = - + J2^A(t), (19a) 
z t=l 

ta = lim ta(w). (19b) 

w— >oo 

The integrated autocorrelation time determines the size of the standard error in measuring 
A according to, 

8 A = (2var(A)r A /iV MC ) 1 / 2 (20) 



with var(A) the variance in A and Nmc the number of Monte Carlo steps. Figure p~5| is a 
plot of ta{w) with A = m, e and / as a function of w for the two-dimensional Ising model. 
Similar behavior was found for the three-dimensional Ising model. The error for ta{w) was 
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estimated by taking the square root of the sum of the variances of r^(t)'s for t < w. For all 
three observables, ta{w) reaches a plateau in a few steps. 

Table | is a summary of the integrated autocorrelation time for the two and three- 
dimensional Ising model at w = 6, where all the r's have saturated but still have relatively 
small errors. The values for r e are compared with results |J for the SW and Wolff single 
cluster algorithm. The energy autocorrelation time is markedly smaller for the IC algorithm 
than the other two cluster algorithm. Furthermore, t £ decreases for larger systems while 
for the other cluster algorithms it increases. r m appears to be independent of system size 
suggesting the possibility that the dynamic exponent for IC dynamics is zero for both two 
and three-dimensional Ising models. 

For the largest systems, 1^(1) is close to —0.5 and, as a result, Tf is very small. Although 
the data is not good enough to draw clear conclusions, it appears that Tf approaches zero 
as a power, L Zf with Zf « — 1 . The anticorrelation in / and e means that these quantities 
can be accurately estimated in a small number of Monte Carlo steps. Indeed, for averaging 
these variables, the IC algorithm is better than performing independent sampling from 
the invariant IC measure. For example, Eq. ^D] and the behavior of var(/), implies that 



6f ~ L- a /^N^ where a « 1.3. 

Table |T| shows results for the integrated autocorrelation time for the 3 and 4-state Potts 
models. We find again that the IC algorithm is much faster than the SW algorithm though 
for q = 3 and 4 there appears to be some critical slowing. Based on Table || we can obtain 
estimates for the dynamic exponet for the magnetization, z m (3) = .28 and z m (4) = .63. 
Note however that these values are less than the Li and Sokal |35j bound for the dynamic 
exponent for the SW algorithm (z > ajv). 

V. SUMMARY 

The invaded cluster method comprises a class of algorithms for sampling equilibrium 
spin systems. Because cluster growth is controlled by a spanning rule rather than the 
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temperature, the method is able to simulate the phase transition point without a priori 
knowledge of the phase transition temperature. The transition temperature is, instead, an 
output of the algorithm. We have demonstrated this numerically for Ising/Potts models in 
two and three dimensions. 

We may also use parameterized stopping rules to explore either the coexistence region 
of discontinous transitions or the critical region near a continuous transition. For these 
rules we specify a quantity such as the energy or susceptibility and learn the corresponding 
temperature. Using the mass rule we have been able to sweep through the coexistence 
region of first-order transitions and to obtain quantities such as the energy of the high and 
low temperature coexisting phases. The behavior of the effective transition temperature with 
the mass parameter apparently yields a very sensitive method to distinguish continuous from 
discontinuous transitions. 

The invaded cluster algorithm is very similar to the Swendsen-Wang algorithm except 
that the occupied bonds are determined by a stopping rule rather than the temperature. 
We argued that this leads to a feedback mechanism that forces the system to the desired 
equilibrium state much faster than is the case for the Swendsen-Wang algorithm. The conse- 
quence is that the algorithm is extremely fast. Measured autocorrelation times are less than 
unity and decrease with system size for the energy and estimated critical temperature. The 
magnetization integrated autocorrelation time is constant for the two and three-dimensional 
Ising models but grows slowly for the 3 and 4-state two-dimensional Potts models. We 
speculate that the invaded cluster algorithm applied to Ising critical points has no critical 
slowing down. For this reason and because there is no need to know the transition tem- 
perature, we believe the IC method will prove to be the most efficient approach for high 
precision measurements of critical properties. 

Although we have tested the algorithm in a number of setting and supplied non-rigorous 
arguments for its validity much work remains to be done in understanding the method and 
putting it on a firm footing. It is important to prove that the IC ensemble is equivalent to 
the usual statistical mechanics ensembles. We would also like to understand the finite-size 
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scaling properties of the IC ensemble since these differs from our naive expectation in some 
cases. 

In this paper we have confined our attention to Potts models however the method is much 
more broadly applicable. In a future paper [[HJ we will show how to use the approach for a 
variety of discrete spin systems such as the Ashkin- Teller model. Similary, the embedding 
approach described by Wolff can be used to generalize the method to 0(n) models. 

This work was supported by NSF Grants DMR-93-11580, DMR-95-0013P and DMS-93- 
02023. 



NOTE ADDED 



After completing this research we received an interesting preprint |36| that describes 
a "fixed cluster" algorithm. This algorithm uses a stopping rule based on the extent I of 
the largest cluster. However, in contrast to the fixed parameter rules used here, I does not 
correspond to a thermodynamic quantity. 
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TABLES 

TABLE I. Integrated autocorrelation times for two- and three-dimensional Ising models for the 
SW, Wolff and IC algorithms. Results for the IC algorithm are measured at time step w = 6. 



d 


L 


T £ ,SW* 


T £ ,Wolff a 


T £ ,IC 




T f,ic 


2 


32 


4.13(4) 


1.80(1) 


.51(3) 


.88(3) 


.19(3) 


2 


64 


4.92(8) 


2.23(3) 


.42(2) 


.78(3) 


.11(2) 


2 


128 


6.00(8) 


2.69(4) 


.42(3) 


.80(3) 


.07(3) 


2 


256 




3.17(8) 


.37(3) 


.81(3) 


.06(3) 


3 


16 


5.6(1) 


1.36(2) 


.35(2) 


.65(2) 


.09(2) 


3 


24 


6.8(1) 


1.50(3) 


.27(2) 


.62(3) 


.07(2) 


3 


32 


7.8(3) 


1.72(4) 


.25(2) 


.65(2) 


.05(2) 


3 


48 


9.9(4) 


1.90(6) 


.19(2) 


.66(3) 


.02(3) 



a Ref. @. 
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TABLE II. Integrated autocorrelation times for two-dimensional, 3 and 4-state Potts models 
for the SW and IC algorithms. Results for r e jc an d t//c f° r IC dynamics are measured at time 
step w = 6 while the time step w m for r m ,ic is shown in the last column. Results for SW dynamics 
are for sizes 128 and 256 rather than 120 and 250. 



q 


L 


T £ ,SW h 


1~e,IC 


T f,IC 


1~m,IC 


W m 


3 


120 


30.3(1.2) 


.73(3) 


.08(3) 


1.40(4) 


6 


3 


250 


39.6(1.7) 


.59(2) 


.06(2) 


1.73(5) 


11 


3 


500 




.52(2) 


.06(3) 


2.08(6) 


15 



4 


120 


115.7(6.1) 


1.23(2) 


.11(3) 


2.97(7) 


16 


4 


250 


232.0(24.6) 


1.10(3) 


.09(2) 


4.61(10) 


27 


4 


500 




0.88(3) 


.05(2) 


7.31(16) 


57 



3 Ref. |5|. 
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FIGURES 

FIG. 1. (/) vs. L~ x l v for the three-dimensional Ising model using different stopping rules. The 
infinite volume estimate of p(T c ) from Ref. |31] is shown on the vertical axis. 



FIG. 2. Double logarithmic plot of {f)t — (f) e vs. L for the three-dimensional Ising model. 

FIG. 3. Double logarithmic plot of the distribution of cluster sizes n(s) for the three-dimensional 
Ising model. The inset shows a double logarithmic plot of the average size of the largest cluster. 

FIG. 4. The mean and the median of / vs. 1/L for the two-dimensional Ising model. The solid 
line shows a linear fit to the median and the exact infinite volume value is shown on the vertical 
axis. The inset shows the standard deviation of the / vs. 1/L. A least squares fit to the form 
co + ciL -1 / 2 + (solid line) suggests that the distribution becomes sharp in the infinite volume 

limit. 

FIG. 5. (e) vs. 1/L for the two-dimensional Ising model. The solid line is a fit to the form 
£o + £iL -1 + £2L~ 2 and the exact infinite volume result is shown on the vertical axis. The inset 
shows var(e)iV vs. L with a linear fit through the data. 

FIG. 6. The mass rule applied to the two-dimensional Ising model. The upper graph shows (/) 
vs. m. In the lower graph the same data is scaled as described in the text. 

FIG. 7. Double logarithmic plots of \ f —p(T c )\ vs. system size L for the two-dimensional g-state 
Potts models. Exact values of p{T c ) are used. 

FIG. 8. Double logarithmic plots of the standard deviation of /, var(/) 1//2 vs. system size L for 
the two-dimensional g-state Potts models. 

FIG. 9. The mean and median of / for the three-dimensional 3-state Potts model. The upper 
graph shows results from the topological stopping rule, the solid line is a linear fit to the median. 
The lower graph shows results from the extension stopping rule. 
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FIG. 10. (/) vs. m for the two-dimensional 10-state Potts model using the mass rule. The 
dashed line marks the exact value, p(T c ). The inset shows the standard deviation of / vs. 1/L for 
m = 0.4 (squares) and m = 0.85 (triangles) . 

FIG. 11. (e) vs. m for the two-dimensional 10-state Potts model using the mass rule for several 
lattice sizes L. The solid line is a fit to the form cq + c\m + C2m 2 for the L = 200 data points with 
0.05 < m < 0.5, the intercept with the m = axis is our estimate for e^, the energy of the high 
temperature phase at the transition . 

FIG. 12. The median of / vs. m for the two-dimensional Potts models using the mass rule with 
q = 2, 4, 5, 6, and 10 and L = 200. The exact value of p(T c ) for each q is shown by a dashed line. 

FIG. 13. / vs. m for the three-dimensional 3-state Potts model. 

FIG. 14. Autocorrelation functions of the magnetization m, energy e and occupation fraction 
/ vs. time step t for the two-dimensional Ising model for size, L = 256. 

FIG. 15. Integrated autocorrelation times for m, e and / vs. integration time, w for the 
two-dimensional Ising model for size, L = 256. 
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